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In Memoriam 

Michael D. Gibbons 
(1960 — 1994 ) 

This report is being published posthumously and is dedicated to the author, 
Michael D. Gibbons. It should be noted that this report was not complete at 
the time of Mike’s death. He had mentioned to several of his co-workers that 
he had gathered all of the information necessary to complete the report, but 
obviously had not had time to process and incorporate all of the data. Mike 
was a meticulous note-taker, and so with very little effort, we were able to find 
and extract the necessary information to complete and publish this valuable, well- 
presented, research. This report stands as a testament to the accuracy, attention 
to detail, and overall thoroughness Mike applied to all facets of his work. 

In the fall of 1994, Mike’s Lockheed and NASA co-workers had resigned 
themselves to the fact that their day-to-day relationship with him was about 
to change drastically, since he had recently accepted a position with Boeing in 
Seattle, Washington. Nobody truly felt Mike was leaving since he would always 
be a phone call, or as he preferred, an e-mail message away. Everyone was 
excited for him and preparing to offer their best-wishes and farewells when he 
tragically died just two days before his scheduled departure. 

Mike dedicated virtually his entire career to aeroelastic analysis and the 
prediction of flutter. He was responsible for the development and verification 
of many of the flutter prediction techniques using computational fluid dynamics 
in practice today. His work in the application of Transonic Small Disturbance 
(TSD) potential flow methods to problems in aeroelasticity is unsurpassed. Due 
to this work, he had established himself as one of the premier computational 
aeroelasticians in the industry. The research discussed in this report breaks ground 
in several different areas. From a numerical standpoint, it is a comprehensive 
study of the impact of varying equation accuracy on unsteady aerodynamic and 
aeroelastic simulations. Data are presented and compared using methods ranging 
from linear theory to TSD to full Euler/Navier-Stokes formulations. Second, and 
more importantly, it provides valuable insight into the role of viscous effects in 
flutter, especially in the vicinity of the “transonic dip.” 

This effort was basically Mike’s first foray into the world of Euler/Navier- 
Stokes aerodynamics and its application to aeroelastic problems. As evidenced by 
this report, it is obvious that he possessed the talent, willingness and dedication 
required to excel in this field, and the industry will miss him dearly. 

To many of us, Mike was a loyal friend who cheered us with his kindness 
and warm smile. He touched us deeply, and along with his family, we grieve 
his untimely death. We will always have those happy memories of working and 
playing with this special friend. His passing has left not only a technical void in 
the Langley community, but an emotional one as well. 








Aeroelastic Calculations Using CFD 
for a Typical Business Jet Model 
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Summary 

Two time-accurate Computational Fluid Dynamics (CFD) codes were used 
to compute several flutter points of a typical business jet wind tunnel model. 
The model consisted of a rigid fuselage with a flexible semispan wing and was 
tested in the Transonic Dynamics Tunnel at NASA Langley Research Center 
where experimental flutter data were obtained from Moo =0.628 to M =0.888. 
The computational results were computed using CFD codes based on the inviscid 
TSD equation (CAP-TSD) and the Euler/Navier-Stokes equations (CFL3D-AE). 
Comparisons are made between analytical results and with experiment where 
appropriate. The results presented here show that the Navier-Stokes method is 
required near the transonic dip due to the strong viscous effects while the TSD 
and Euler methods used here provide good results at the lower Mach numbers. 

Nomenclature 
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C r 


Ct 
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f, Hz. 

fi 
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Moo 

<li 

Qoo 

Q 
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local chord 
root chord (1.8942 ft) 
tip chord (0.5433 ft) 
damping matrix 
sectional lift coefficient 
pressure coefficient 
frequency 

mode shapes (i=l, 6) 
stiffness matrix 
mass of wing (0.56522 Slugs) 
mass matrix 

freestream Mach number 
generalized displacements (i=l, 6) 
freestream dynamic pressure 
generalized forces 
semi-span of wing (4.0392 ft) 
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t nondimensional time 

Uf flutter velocity 

V volume of truncated cone that encloses wing (5.1945 ft 3 ) 

Vf flutter speed index 

x,y,z nondimensional cartesian coordinates 

7 ratio of specific heats 

fi mass ratio 

pi flutter density 

<f> small disturbance potential 

ui flutter frequency, rad./sec. 

uji modal frequencies (i=l, 6), radVsec. 


Introduction 


Several years ago a significant effort was devoted to the development of a 
Computational Fluid Dynamics (CFD) method that solves the inviscid Transonic 
Small Disturbance (TSD) equation for use in aeroelastic analysis. This work was 
undertaken in order to provide an alternative to the widely used linear methods, 
such as doublet lattice, and to improve upon previous TSD codes. The code 
developed to solve the TSD equation is called CAP-TSD 11, 21 and allows transonic 
effects to be included in flutter calculations. Subsequent to the CAP-TSD effort, 
a version of CFL3D^ 3, which solves the Euler/Navier-Stokes equations was 
modified to allow dynamic aeroelastic calculations.^^ This version of the code is 
known as CFL3D-AE, and it was used to compute all of the Euler/Navier-Stokes 
calculations presented in this report. The advantage of using the Euler/Navier- 
Stokes equations is the improved accuracy in modeling the flow physics. Both 
CAP-TSD and CFL3D were developed at NASA Langley Research Center. 

The purpose of this report is to further assess the above CFD methods for 
transonic flutter applications. Previously CAP-TSD has been used to compute the 
flutter boundary of various wings, including thin highly swept delta wings J 6, ^ 
the AGARD 445.6 wingj 8 ^ and an Active Flexible WingJ 9 ^ CFL3D has been used 
to solve the Euler^ and the Navier-Stokes^ 1 ^ equations for the AGARD 445.6 
wing. This report presents results for a typical business jet configuration 112 ^ that 
was tested in the Transonic Dynamics Tunnel at NASA Langley Research Center. 
The results include spatial and temporal convergence studies, surface pressure 
coefficient comparisons for rigid and statically deformed cases, and a discussion 
of the aeroelastic results using the TSD, Euler, and Navier-Stokes methods. 
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Numerical Method 


This section gives a brief description of the methods used to numerically 
solve the TSD, Euler/Navier-Stokes, and aeroelastic equations as contained in the 
CAP-TSD and CFL3D computer codes, respectively. 

The CAP-TSD code solves the unsteady transonic small disturbance equa- 
tion using an implicit time-accurate approximate factorization algorithm with the 
optional use of internal subiterations to minimize linearization and factorization 
errors. The unsteady aerodynamics are simultaneously integrated in time along 
with the structural equations of motion, which allows a time history of the struc- 
tural response to some initial conditions to be studied. The structure is modeled as 
a series of orthogonal mode shapes weighted by time varying coefficients known 
as the generalized displacements. Since the modal deflections in the streamwise 
and spanwise directions are small in comparison to the vertical modal displace- 
ments, they are neglected and the position of the wing at any point in time can 
be given by 


modes 

z(x,y,t)= 52 qi(t) f i( x >y) 

i 


where qj is the time varying generalized displacements and £ represents the 
vertical components of the mode shapes. This allows the structural equations 
of motion to be written in generalized coordinates 


[M]{q} + [C]{q} + [K]{q} = {Q} 


where [M], [C], and [K] are the generalized mass, damping and stiffness matrices 
and {Q} is the generalized aerodynamic forces. The unsteady aerodynamics 
can be computed using different forms of the TSD equation by choosing different 
coefficients. Two different forms of the TSD equation were used here by choosing 
either the linear equation coefficients or the so-called AMES coefficients. The 
linear results were computed using the form 


j t (-MUt - 2 Ml*,) 
+ ±((1 - Ml )*.) 
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The nonlinear results were computed using the TSD equation written using the 
AMES coefficients given by 


If (-■*£*- 2A4^) 

+ £((i - *&)*« - i ( 7 + + “(7 - wifi ) 

+ ^(« = 0 

When the linear equation was used, the wing was modeled as a flat plate in order 
to produce results similar to other linear methods such as doublet lattice. When the 
nonlinear equation was used, the wing was modeled using the appropriate airfoil 
shapes so that nonlinear effects (such as moving shock waves) could be modeled. 
All CAP-TSD calculations included the effects of shock generated entropy and 
vorticityJ 13 ^ The linear pressure coefficient equation was used to compute the 
surface pressure values. Version 1.3 of CAP-TSD was used in all calculations. 

The CFL3D code solves either the Euler or Navier-Stokes equations using an 
implicit time-accurate approximate factorization algorithm involving either flux 
vector splitting or flux difference splitting. For the calculations presented here, 
flux-vector splitting was used. For the Navier-Stokes calculations, turbulence was 
modeled using the Baldwin-Lomax model. Coupling of the structural equations of 
motion with the aerodynamics was done similar to that implemented in CAP-TSD. 
As with CAP-TSD, only the vertical component of the mode shape is used. 


Aeroelastic Calculations for a Typical Business Jet 

This section discusses the numerical results obtained for a business jet con- 
figuration. Details of the geometry and the CFD models used in the calculations 
are included. Every attempt was made to verify the spatial and temporal con- 
vergence of the solutions. Results are presented which show the convergence of 
selected cases. Sensitivity of the aeroelastic solutions to the initial conditions are 
also presented. Surface pressure comparisons are made for the static rigid and 
static aeroelastic cases. Finally, dynamic aeroelastic results are presented using 
different equation levels, and comparisons are made with experimental data. 


Model Geometry 

The business jet configuration consisted of a semispan wing/fuselage con- 
figuration which was tested in the Transonic Dynamics Tunnel (TDT) at NASA 
Langley Research Center. Figure 1 shows an illustration of this model mounted on 
the wall of the TDT. The model was tested in air for a range of freestream Mach 
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numbers to obtain flutter points. The model, however, was not instrumented to 
measure surface pressure data. Experimental flutter data were obtained for Mach 
numbers ranging from 0.628 to 0.888 [14] . The wing root angle-of-attack was var- 
ied during the test to minimize loading with the maximum angle needed for this 
purpose being 0.2 degrees. The wing has a taper ratio of 0.29 and a midchord 
sweep of 23 degrees. The airfoil thickness varies from 13 percent at the symmetry 
plane (for the extended wing configuration) to 8.5 percent at the wing tip. The 
fuselage has a circular cross-section with a conical aft end. 

The CAP-TSD calculations for the wing/fuselage configuration involved a 
100x50x80 point computational grid with 45 points along the chord of the 
wing and 25 points along the span. Figure 2 shows a perspective view of the 
wing/fuselage surface grid. Figure 3 shows the four airfoil sections used to define 
the wing surface. Linear interpolation is used to define the airfoil sections at 
other spanwise locations. Calculations were also performed for a wing alone 
configuration to study the effect of the fuselage on the results. Calculations for 
the wing alone used the same grid as the wing/fuselage configuration but with the 
wing extended to the symmetry plane. 

Figure 4 shows the mode shapes for the wing, and Table 1 gives the general- 
ized masses and frequencies used in the aeroelastic calculations. The mode shapes 
and generalized masses were computed from a finite element structural model us- 
ing NASTRAN. The frequencies are experimental values measured during testing 
of the model in the TDT. As with most models of this type, a yaw mode is present 
which includes very small vertical displacements. In all subsequent analyses, this 
yaw mode is neglected and is not shown in the above figure. 

The CFL3D calculations for the extended wing were made primarily with two 
grids as shown in Figures 5 and 6. The Euler mesh contained 153 points around 
the airfoil and wake, 57 points out the span with some clustering around the wing 
tip, and 33 points normal to the wing surface. There were 113 points around each 
airfoil section and 41 points on the wing spanwise. The Navier-Stokes mesh was 
vir tuall y identical except there were 51 points normal to the wing. Both meshes 
used a C-H topology. No attempt was made to model the fuselage since inviscid 
TSD results showed that it has little effect on the flutter solution as discussed 
later. Viscous interactions between the wing and fuselage could impact flutter 
results, but detailed analyses of these types of interactions are beyond the scope 
of this effort. 

Spatial and Temporal Convergence 

Prior to making numerous aeroelastic calculations, a convergence study was 
done with CAP-TSD and the CFL3D-Euler results. The effect of time step 
on the CAP-TSD solutions is shown using the linear and nonlinear equations, 
and the effect of grid density is shown for a case using the nonlinear equation. 
Spatial and temporal convergence of the CFL3D-Euler results are shown. No 
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convergence study was done using the CFL3D-Navier-Stokes version due to the 
large computational time required. By making these calculations, a grid density 
and time step were chosen which produced solutions of acceptable accuracy and 
yet require the minimum amount of computer time. 

CAP-TSD 

The temporal convergence study was performed at Moo=0-6 with CAP-TSD 
using the nonlinear equation and at Moo=0.9 using the linear and nonlinear 
equation. Table 2 shows the effect of time step on the predicted flutter dynamic 
pressure at Moo=0.9 using the linear TSD equation. The table shows the flutter 
dynamic pressure and frequency computed using a time step based on 165, 330, 
and 660 steps per cycle of motion in the third modal frequency (31.6 Hz). 
Increasing the number of time steps from 165 to 330 results in roughly a 7 
percent reduction in predicted flutter dynamic pressure whereas increasing from 
330 to 660 only produces about a 1 percent reduction in flutter dynamic pressure. 
Figure 7 shows the first generalized displacement computed at Moo=0.6 using 
the exact same input files except for the time step. The figure shows the large 
error when the solution is not converged temporally. Using a time step of 0.1 
results in an unstable structural response whereas a value of 0.0339 produces a 
stable structural response. Reducing the time step further to 0.0166 resulted in no 
significant change from the response at 0.0339. The nonlinear case at M<x>=0.9 
is shown in Figure 8. The transients shown correspond to 330 and 660 steps 
per cycle of motion. Figure 8 indicates that at Moo=0.9 a temporally converged 
solution was obtained when using 330 steps/cycle. The time step used in the 
remaining unsteady calculations was based on 330 steps per cycle of motion 
since the increased accuracy at 660 was deemed minimal. 

A spatial convergence study for a nonlinear case at Moo=0.9 and qoo=21.6 
psf is shown in Figure 9 using two different grid densities. The first grid size 
was 100x50x80 and is used in all the other results computed in this paper; the 
second grid size was 140x70x90. The second grid has 60 points along the chord 
and 35 points along the span with clustering near the wing tip. Figure 9 shows 
that the two solutions are nearly identical. This good comparison suggests that 
the 100x50x80 grid produces solutions of acceptable spatial accuracy. 

CFL3D-Euler 

Three different grid densities were used to study spatial convergence of the 
CFL3D-Euler solution at Moo=0.9 and zero degrees angle of attack. The grids 
used were constructed using a C-H topology. These grids contained 153x43x33, 
153x57x33, and 263x57x33 grid points. The index refers to the number of points 
around the airfoil chord, spanwise, and normal to the surface, respectively. Figure 
10 shows the comparisons between the three solutions. Each grid gives nearly 
the same solution with the spanwise index having the greatest influence on the 
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location of the upper surface shock near the wing tip, and the chordwise index 
affecting the resolution of a small lower surface shock near the tip. Based on 
these calculations, the 153x57x33 grid was selected as producing a solution with 
an acceptable level of spatial accuracy for the aeroelastic calculations presented 
in this report. 

The CFL3D-Euler calculations used a time step of 0.05 for all the unsteady 
results presented. To determine how sensitive the solutions were to time step, a 
calculation was done with a time step of 0.01 for an aeroelastic case at M 00=0.9 
and qoo=36.0 psf. Figure 11 shows a plot of the first generalized displacement 
for the two different time steps. Since the two solutions are nearly identical, the 
larger time step of 0.05 was used for further dynamic aeroelastic calculations. 

Initial Conditions 

Each dynamic aeroelastic calculation is computed by using the converged 
static aeroelastic solution with some form of initial condition prescribed for the 
generalized displacements or velocities of each structural mode. Typically, the 
initial condition on the generalized displacements of each mode is taken to be 
zero which should help to reduce numerical transients which otherwise might be 
created by the wing instantaneously being displaced. To initiate the motion of 
the wing, the generalized velocities of each mode have been set to 1.0 in all 
the dynamic aeroelastic calculations presented unless noted otherwise. By setting 
all of the generalized velocities to 1.0 (ie, qi = 1.0, i = 1,6), each mode actively 
participates in the structural response. This allows for accurate values of damping 
and frequency to be extracted from each generalized displacement. 

The choice of using an initial generalized velocity of 1.0 is somewhat arbitrary 
and its effect on the structural response depends on how the mode shapes were 
scaled. To show the effect that different values have on the solution, values 
of 0.1, 1.0, and 10.0 were used as initial conditions at Moo=0.85, qoo=64.8 
psf for the wing/fuselage configuration. These calculations were done using the 
CAP-TSD code. Figure 12 shows the structural response of the first generalized 
displacement. Each response was normalized by the initial condition amplitude 
for direct comparison. Figure 12 indicates that the responses are quite similar, 
and when damping and frequency information is extracted from each curve, the 
frequency values are nearly identical, as shown in Table 3, while the damping 
values range within five percent of each other. The amplitude of motion at the 
wing tip leading edge caused by the different initial conditions varied from a 
maximum of +/- 0.015 degrees to +/- 1.5 degrees for the initial conditions 0.1 
and 10.0, respectively. These calculations indicate the aeroelastic response to 
a set of initial conditions behaves in a linear fashion with the initial conditions 
even though the mean flow field contains nonlinear features. Similar results were 
obtained at a higher Mach number, Moo=0.9, for the wing alone case as shown 
in Figure 13 and Table 4. 
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Methodology to Locate Flutter Crossing 

To compute the point at which flutter occurs for a given Mach number, several 
executions of a code are required. First, a static aeroelastic solution is computed 
using a value of the dynamic pressure that is assumed to be near that which pro- 
duces flutter (neutral stability). Next, a dynamic aeroelastic solution is computed 
by restarting the calculation from the converged static aeroelastic solution with 
some initial condition on the vertical velocity of the wing. The generalized dis- 
placements from the dynamic case are fit with a series of exponentially damped 
sine waves^ 15 ^ allowing the stability of the system to be determined. If the aeroe- 
lastic solution was stable, a higher value for the dynamic pressure is chosen and 
the static and dynamic aeroelastic solutions are again computed. Once two or 
more transients have been computed, the flutter dynamic pressure value is deter- 
mined by linear extrapolation using the damping information from the curve fitting 
process. Linear extrapolation or interpolation gives an estimate of the conditions 
that produce flutter. However, further refinement can be obtained by additional 
aeroelastic calculations if improved accuracy is desired. 

An alternative method to calculate the flutter crossing is to vary free stream ve- 
locity while holding density constant (usually to the experimental flutter density). 
This approach allows the mass ratio of the calculation to match the experimental 
value. Since this method generally requires that the flutter solution be known be- 
fore it is computed, the method of varying density is generally used. An additional 
reason to vary density rather than velocity is that the flutter dynamic pressure is 
generally more sensitive to changes in density than velocity thus making it easier 
to determine the point of neutral stability. 

Static Rigid Results 

In this section, comparisons of the surface pressure coefficients computed 
with CAP-TSD and CFL3D are shown for the extended wing at zero degrees 
angle of attack. Due to the lack of experimentally measured surface pressures, 
comparisons are made between different codes. The most direct comparison to 
make with CAP-TSD is the CFL3D-Euler results since both methods assume the 
flow to be in viscid. Comparisons between CFL3D-Euler and CFL3D-Navier- 
Stokes generated surface pressure distributions are shown to indicate the effect 
that viscosity has on the steady solution. These calculations were made for the 
extended wing configuration only since the fuselage was not modeled in the 
CFL3D calculations. 

Comparisons between CAP-TSD and CFL3D-Euler surface pressure coeffi- 
cients are shown for Moo=0 7 and Moo=0.9 at zero degrees angle of attack. Figures 
14 and 15 show that CAP-TSD compares favorably with the CFL3D-Euler re- 
sults. At Moo=0.7, both methods are in excellent agreement for this subsonic case. 
When the Mach number is increased to 0.9, a shock wave develops along the up- 
per surface with a weaker shock wave on the lower surface. Comparisons between 
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the two methods indicate a slight difference in the surface pressure magnitude at 
some span stations. On the lower surface near midspan, CAP-TSD predicts the 
shock is located further upstream than the Euler results. Part of this difference 
may be due to the CAP-TSD grid which had a finer spacing in the spanwise direc- 
tion near the wing root than the Euler grid did. Although not shown, additional 
CAP-TSD calculations with a finer spacing at the wing trailing edge, resulted in 
excellent comparisons of the surface pressure coefficients near the wing trailing 
edge but the refinement had little effect on the overall surface pressures. 

Comparisons are also shown between the sectional lift coefficients for the 
above two cases in Figure 16. As expected, the comparison between CAP-TSD 
and CFL3D is fairly good at Moo=0.7 but shows some significant differences at 
Mqo =0.9. These differences are due to the slight disagreement in the surface 
pressures shown in Figure 15. 

To show the effect that viscosity has on the position and strength of the shock 
waves, a comparison was made between the CFL3D-Euler and CFL3D-Navier- 
Stokes results at Moo=0.9, zero degrees angle of attack, and Re = 1.107xl0 6 . This 
Reynolds number is typical of the value obtained during flutter tests this model 
in the TDT. Figure 17 shows how the shock on the upper surface is significantly 
weakened and shifted forward. The rapid recompression at the trailing edge is 
also weakened due to the presence of the boundary layer. Aeroelastically, the 
most important features are the forward shift in the shock resulting in a forward 
shift in the aerodynamic center and the lower loads experienced on the wing due 
to the weakening of the shock. 

Static Aeroelastic Results 

Prior to making a dynamic aeroelastic calculation, the wing was allowed to 
deform due to the static air loads. To obtain a converged static aeroelastic solution, 
six normal modes were used. A structural damping value of 0.99 (nearly critically 
damped) was used for each mode while the solution was marched in time. Surface 
pressure distributions are shown, and comparisons are made between CAP-TSD 
and CFL3D-Euler results. The aerodynamic center and its relationship to the first 
torsion node line, based on CAP-TSD calculations, are also plotted. 

Figure 18 shows a three dimensional view of the upper and lower surface 
pressure coefficients computed with CAP-TSD for the wing/fuselage configura- 
tion. Results are shown at Moo=0.8, 0.85, 0.9, and 0.95 for the converged static 
aeroelastic deflections near the flutter dynamic pressure. These figures show the 
development of strong shocks on the surface of the wing. At Moo=0.95, the upper 
and lower surface shock has reached the trailing edge with the recompression oc- 
curring slightly downstream of the wing trailing edge. The occurrence of a shock 
near the wing root trailing edge is due to the aft shape of the fuselage which 
results in a rapid expansion followed by a shock. As Figure 2 shows, the fuselage 
becomes conical at approximately eight tenths root chord. Figure 18 shows that 
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the shock completely disappears when the aft end of the fuselage was replaced 
by a cylinder extending to the downstream boundary. In reality, the development 
of such strong shocks may not occur due to the effect of viscosity which tends 
to weaken shocks and move them forward. 

Comparisons between the CAP-TSD and CFL3D-Euler static aeroelastic so- 
lutions at Moo =0.7 and Moo=0.9 are shown in Figures 19 and 20. At Moo =0.7 
and qoo=l 15.2 psf, comparisons between the CAP-TSD and CFL3D-Euler results 
are quite good. At this Mach number a strong suction peak exists along the 
lower surface leading edge. This is due to the negative effective angle of attack 
produced by the twist of the deformed wing. 

At Moo =0.9 and qoo=36.0 psf, comparisons between CAP-TSD and CFL3D- 
Euler pressures (Figure 20) are good. However, CAP-TSD predicts the lower 
surface shock wave to be located further forward near the wing tip than predicted 
by CFL3D. Part of the difference in the solutions is due to the statically deformed 
shape. The loading predicted with CAP-TSD results in the wing tip leading 
edge being deflected 1.81 degrees nose down whereas CFL3D-Euler predicts a 
deflection of 1.74 degrees. Since the loads on the wing were slightly different for 
the rigid case, the resulting static aeroelastic shape must be different. 

A comparison of the sectional lift at Moo=0.7 and Moo=0.9 is shown in Figure 
21. At Mqo =0.7, the comparison between CAP-TSD and CFL3D sectional lift is 
fairly good with some discrepancy near the wing tip. At the higher Mach number 
of Moo=0.9, CAP-TSD tends to predict a higher loading at the wing tip and a 
lower loading at wing root again due in part to the difference in the statically 
deformed shape. 

Figure 22 shows how the aerodynamic center, computed with CAP-TSD near 
the flutter dynamic pressure, varies as the Mach number increases from Moo=0.7 
to Mqo =0-9. At Moc=0.7, the aerodynamic center is located well forward of the 
first torsion node line near the quarter chord line as expected. This is one reason 
why subsonically a higher dynamic pressure is required to produce flutter. As 
the Mach number increases, the aerodynamic center moves aft and at Mqo =0.9 
lies close to the first torsion node line. This reduces the aerodynamic damping 
of the wing. With less damping the dynamic pressure which will produce flutter 
is expected to decrease. Although not shown, calculations where the dynamic 
pressure was varied at a constant Mach number indicated little change in the 
location of the aerodynamic center. 

Dynamic Aeroelastic Results 

The flow conditions which produce flutter that are presented in this section 
were computed using TSD, Euler and the Navier-Stokes equations. Calculations 
for the wing/fuselage configuration and the extended wing configuration were 
made with CAP-TSD. Since the presence of the fuselage appeared to have little 
effect on the aeroelastic solutions, the Euler and Navier-Stokes calculations were 
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made only for the extended wing. Calculations are presented to show the effect 
of the fuselage, mass ratio, and equation level. 

For each calculation the results were tabulated and plotted. To eliminate any 
possible confusion on how certain quantities were calculated, the expressions used 
are given here. The mass ratio was computed using 


A* = 


m 

Wf 


where 

V = — Jra(c? + c t c T + c 2 r ) 

and ‘m’ is the mass of the wing, ‘s’ is the semispan of the wing, ‘c<’ is the tip 
chord, and ‘c r ’ is the root chord at the wing fuselage intersection. The subscript on 
density refers to the value at flutter. The reduced frequency was computed using 

k= 2U7 


where u>{ is the flutter oscillatory frequency, and Uf is the flutter free stream 
speed. The flutter speed index was computed using 


V f 


2U f 

C t ^3 y/ji 


where ujz is the third modal frequency corresponding to first torsion. 


CAP-TSD results 

Several flutter points were computed with CAP-TSD using the linear and 
nonlinear equations. For both sets of results, the flutter points were computed by 
holding velocity constant and varying density unless stated otherwise. No struc- 
tural damping was used in these calculations. All calculations were performed 
with the wing root set at zero degrees angle of attack. Each dynamic aeroelastic 
calculation was restarted from a converged static aeroelastic solution. The ma- 
jority of the unsteady calculations were computed using 330 steps per cycle of 
the third modal frequency. Typical CPU execution times were 3.1 micro seconds 
per grid point per time step on a CRAY 2. Flutter data for the business jet wing 
are tabulated in Tables 5 through 13. Experimental data obtained in the NASA 
Langley Research Center Transonic Dynamics Tunnel are presented in Table 5 
for comparison with the following computed results. 

Figure 23 (Table 6) shows that the CAP-TSD linear results for the 
wing/fuselage configuration at low Mach number are in excellent agreement with 
the experimental results. The good agreement between CAP-TSD linear results 
and experiment at the lower Mach numbers should be expected since nonlin- 
ear aerodynamic effects are insignificant whereas linear theory is thought to give 
reasonable answers at higher Mach numbers due to the counteracting effects of 
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thickness and viscosity. Thickness has a destabilizing effect at transonic speeds 
for inviscid flows whereas viscosity tends to be stabilizing. A comparison between 
the CAP-TSD results computed with the linear (Table 6) and nonlinear (Table 7) 
equation in Figure 23 show that thickness has a significant effect above Moo =0-7. 
The large differences between the results using the linear and nonlinear equations 
are due to the presence of strong shocks in the flow field. 

For the results computed with the nonlinear equation, calculations at Moo =0.9 
indicated a very low dynamic pressure may have been needed to obtain a flutter 
crossing. At this Mach number all solutions computed by varying density down 
to a dynamic pressure of 7.2 psf were unstable with low damping. Since such 
an unrealistically low dynamic pressure may have been needed to obtain a flutter 
crossing no point was plotted for this case. 

The effect of the fuselage on the flutter results is shown by two calculations 
in which the fuselage was removed and the wing extended to the symmetry plane. 
Results of this analysis are tabulated in Table 8. As shown in Figure 24, the effect 
of the fuselage is minimal. This is expected since the generalized aerodynamic 
forces are weighted by the mode shapes, and near the root where the interference 
effects of the fuselage are greatest, the modal displacements are nearly zero thus 
resulting in a contribution to the generalized aerodynamic forces which is far 
smaller than the contributions from the loading near the wing tip where the modal 
displacements are greatest. 

Figure 25 (Table 9) shows the effect of matching the mass ratio to the 
experimental value. These calculations are done by holding density constant 
at a linear extrapolation of the experimental density and varying both dynamic 
pressure and velocity. At Moo=0.8, there is no significant difference between 
the calculations which match velocity or density. However, at Moo=0.9, a solid 
flutter point was obtained with CAP-TSD while the velocity was being varied. 
This is in contrast to the unstable solutions computed by varying density while 
holding velocity constant. Based on the improved solutions computed by varying 
velocity, one should not conclude that this provides more accurate answers than 
varying density. Varying velocity effectively changes the phase lag, which has 
been pointed out by Zwann [16 ^ as being a primary cause of the transonic dip. Thus, 
if the theoretical flutter velocity is significantly different from the experimental 
velocity, the flow physics may not be the same as those in the experimental data 
making direct comparisons meaningless. This implies the velocity should be held 
fixed at the experimental value; however, a somewhat different argument might 
be made for varying density which affects wing loading. 

Figure 26 provides further insight into the inviscid flutter by showing the 
modal interactions for the previous CAP-TSD results computed with the nonlinear 
equation for the wing/fuselage configuration in root locus form. These figures 
show frequency versus damping for each mode with stable roots on the left and 
unstable roots on the right of each of the figures. In Figures 26a, 26b, and 26c, the 
mode which first becomes unstable is first bending. These results are very similar 
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to those obtained by Mohr, Batina, and Yang^ 17 ^ who show a study of the flutter 
mechanism for a wing similar to the one studied here. In Reference [17] however, 
the mass ratio was held constant, whereas in the present study, the mass ratio was 
allowed to vary with increasing Mach number. Calculations with CAP-TSD at 
lower mass ratios for the Moo =0.9 case (Table 10) showed that as the mass ratio 
decreases, the dynamic pressure which produces flutter increases, and the type of 
mode shift shown in Reference [17] occurs; that is, the first torsion mode appears 
to be the primary flutter mechanism. With the higher mass ratios encountered for 
the CAP-TSD calculations varying density, there does not appear to be a shift 
in the primary flutter mechanism, which remains first bending into supersonic 
flow. Reference [17] attributes the mode shift to the upper surface shock passing 
across the first torsion node line. Similar results shown here indicate that as 
the shock moves aft, there is a reduction in aerodynamic damping due to the 
aft movement of the aerodynamic center. As the Mach number increases from 
Moo =0.7 to Moo =0.85, the largest stable value of damping gradually decreases 
with the frequency decreasing for zero damping. The decrease in frequency is 
also shown in Figure 23b. Further calculations have shown that the third mode 
eventually goes unstable, but at much higher dynamic pressures. At Moo =0.9, both 
the first bending mode and first torsion go unstable and appear to be unstable even 
at dynamic pressures near zero (large mass ratios). 

CFL3D results 

The following section presents aeroelastic results which were computed using 
the Euler and Navier-Stokes equations as solved by CFL3D. These results were 
computed using the extended wing configuration with no structural damping and 
with the wing root angle-of-attack at zero degrees. Both the Euler and Navier- 
Stokes solutions were computed by holding velocity constant and varying density 
unless stated otherwise. As with the TSD calculations, the dynamic aeroelastic 
calculations were restarted from a converged static aeroelastic solution. The 
Euler calculations were made with a 153x57x33 grid while the Navier-Stokes 
calculations used a 153x57x51 grid. Typical execution times were approximately 

46.0 micro seconds per grid point per time step for the Euler calculations and 

50.0 micro seconds per grid point per time step for the Navier-Stokes calculations 
on a CRAY 2. 

Figure 27 shows a comparison between Euler (Table 11) and TSD (Table 8) 
results. The comparison shows that both methods give nearly the same result 
at Moo=0-7 and Moo=0.85. At slightly higher Mach numbers of 0.92 and 0.94, 
the Euler calculations predict unstable responses down to 7.2 psf similar to the 
results at Moo=0-9 computed with CAP-TSD. Apparently, inviscid solutions for 
this configuration do not contain sufficient aerodynamic damping and result in a 
free response as the Mach number approaches values around 0.92. Both Euler 
and TSD calculations performed near this Mach number exhibit very low flutter 
dynamic pressures, which are not displayed on this plot. As the Mach number is 
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increased to Moo=1.05, a clear boundary between stable and unstable aeroelastic 
transients allowed a flutter point to be computed. An examination of the frequency 
at flutter indicates a mode shift has occurred which is apparently due to second 
bending. 

Figure 28 shows the effects of mass ratio on the computed flutter points. 
Since varying velocity while holding density constant had a significant effect on 
the CAP-TSD solution at Moo=0.9, this method was used to compute the flutter 
crossing with CFL3D-Euler (Table 12). As with the CAP-TSD mass ratio results 
shown in Figure 25, there is little difference in the method used to compute 
the flutter point below Moo=0-85. At Moo=0.9 the method of varying velocity 
predicts a slightly higher flutter speed index. 

The effects of viscosity on the flutter solution are shown in Figure 29. These 
calculations were computed by holding velocity constant and varying density using 
the Navier-Stokes equations. Figure 29 compares the Navier-Stokes solutions with 
the Euler results. Comparisons show that the addition of viscosity results in a small 
increase in the flutter speed at Moo =0-85 while at Moo =0.9 there is a significant 
increase in the predicted flutter speed index. As the Mach number is increased 
further to Moo =0.92, the predicted boundary begins to turn upwards indicating 
that the transonic dip occurs near Moo=0.9. The significant difference between 
the Euler solution at Moo =0.9, which is effectively at an infinite Reynolds number, 
and the Navier-Stokes solution indicates the importance that viscosity plays for 
this configuration. Since the Reynolds number is approximately one million 
for both the Navier-Stokes calculations and experimental data near Moo=0.9, 
it is interesting to note that the experimental data may be nonconservative in 
comparison to flight tests since they usually occur at Reynolds numbers greater 
than one million. The Navier-Stokes results presented in Figure 29 are tabulated 
in Table 13. 


Conclusions 

Detailed aeroelastic calculations have been presented for a business jet config- 
uration computed using the TSD equation, Euler equations, and the Navier-Stokes 
equations. The TSD calculations were computed with CAP-TSD while the Euler 
and Navier-Stokes were computed with CFL3D. 

Comparisons with experiment show that CAP-TSD does a fairly good job 
in predicting the flutter boundary for the business jet model at the lower Mach 
numbers with the results becoming conservative for the higher Mach numbers. 
The reason for the discrepancy appears to be due to the inviscid nature of the 
calculations. Without viscosity, shocks are typically much stronger and further 
aft resulting in greater loads at a given value of dynamic pressure. Calculation of 
the flutter points by varying density or velocity showed no significant differences 
although varying velocity at Mco^.9 did allow for a flutter point to be obtained. 
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Calculations with the Euler equations showed little difference from the CAP- 
TSD solutions; both methods predicted a rapidly dropping flutter boundary as 
Moo=0.9 was approached. Apparently, the assumptions made with TSD, such as 
small perturbations and irrotational flow, are adequate for these cases. Calcu- 
lation of the flutter points by varying density or velocity showed no significant 
differences. 

The results computed with the Navier-Stokes equations showed a significant 
improvement in the comparisons between the computed flutter points and experi- 
ment. The remaining differences between experiment and computed values might 
be further improved by using a turbulence model which models separation such as 
the Johnson-King turbulence model. Some improvement might also be obtained 
by using a finer mesh. Other factors affecting the accuracy of the computed 
flutter points include the accuracy of the structural mode shapes. The structure 
is generally modeled using a subset of orthogonal mode shapes and an accurate 
correlation of the computed mode shapes with experimentally measured values 
may be necessary; if, for instance, the torsion node line near the wing tip were 
inaccurately predicted, it could result in a completely different pitching moment. 
Additionally, the structural damping of the wing should be measured and included 
since the transonic dip region tends to be a region of low aerodynamic damping. 

The goal of this report was to further study the use of CFD for predicting 
flutter. The results presented here indicate that good agreement with experimental 
data can be obtained with the Navier-Stokes code near the transonic dip. At 
the lower Mach numbers, the TSD and Euler results show good agreement with 
experimental data. Comparisons between the various equations levels showed the 
importance of modeling viscosity near the transonic dip and how it had little effect 
on the results at the lower Mach numbers. Since different configurations may be 
more or less sensitive to viscous effects, then obtaining good comparisons with 
experimental data at a minimum cost (CPU time) requires selecting the appropriate 
method for the flow conditions. 
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Table 1 Wing mode shape information. 


Mode Shape 

Generalized Mass 

/]bf sec 2 \ 

^ inch ' 

Frequencies (Hz.) 

1 

0.0040196 

4.30 

2 

0.0019749 

14.50 

3 

0.0049215 

31.60 

4 

0.0012095 

33.00 

5 

0.0018510 

58.20 

6 

0.00060438 

62.66 


Table 2 Effects of number of time steps per cycle of motion on the dynamic pressure and 
frequency at flutter computed with CAP-TSD using the linear equation at M oo =0.9. 


steps/cycle 

165 

330 

660 ; 

dynamic pressure (psf) 

91.38 

85.23 

84.25 

frequency (Hz) 

9.37 

8.82 

8.74 


Table 3 Effects of generalized velocity initial conditions on the damping 
and frequency content of the first generalized displacement computed with 
CAP-TSD using the nonlinear equation at Moo=0-85 and qoo=64.8 psf. 



4i = o-i 

4i = 1.0 

g, = 10.0 

damping 

-0.00477 

-0.00466 

-0.00454 

frequency 

7.972 

7.971 

7.975 


Table 4 Effects of generalized velocity initial conditions on the damping 
and frequency content of the first generalized displacement computed with 
CAP-TSD using the nonlinear equation at Moo =0 90 and qoo=21 .6 psf. 



4i = i.o 

• 55 . 

II 

O 

b 

damping 

-0.0117 

-0.0118 

frequency 

5.603 

5.602 
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Table 5 Experimental flutter data for the business jet wing obtained in 
air in the NASA Langley Research Center Transonic Dynamics Tunnel. 


Mach 

No. 


Dynamic 

Pressure 

(psf) 

Density 

(Slug/ft 3 ) 

10" 4 

Freq. 

(Hz) 

Reduced 

Freq. 

Mass 

Ratio 

0.628 

708.1 

148.3 

5.91 

13.1 

0.110 

184.1 

0.698 

781.0 

140.6 ! 

4.61 

12.5 

0.095 

236.0 

0.781 

863.2 

125.3 

3.36 

11.3 

0.078 

323.8 

0.844 

929.2 

109.5 

2.54 

10.4 

0.067 

428.4 

0.888 

971.2 

106.8 

2.26 

9.3 

0.057 

481.5 


Table 6 CAP-TSD flutter results computed using the linear equation for the wing/fuselage configuration. 


Mach 

No. 

Velocity 

(ft/sec) 

Dynamic 

Pressure 

(psf) 

Density 

(Slug/ft 3 ) 

1(T 4 

Freq. 

(Hz) 

Reduced 

Freq. 

Mass 

Ratio 

Flutter 

Speed 

Index 

0.600 

669.6 

146.04 

6.51435 

12.96 

0.1152 

167.03 

0.2755 

0.700 

781.2 

134.57 

4.41015 

12.08 

0.0920 

246.73 

0.2645 

0.800 

880.0 

117.23 

3.02763 

10.83 

0.0732 

359.40 

0.2468 

0.850 

935.0 

104.03 

2.37990 

9.98 

0.0635 

457.21 

0.2325 

0.900 

990.0 

85.23 

1.73921 

8.82 

0.0530 

625.63 

0.2105 

0.980 

1078.0 

31.80 

0.54729 

5.92 

0.0327 

1988.2 

0.1286 

1.050 

1155.0 

91.33 

1.36924 

9.21 

0.0475 

794.7 

0.2176 


Table 7 CAP-TSD flutter results computed using the nonlinear equation for the wing/fuselage configuration. 


Mach 

No. 

RtfHj 

Dynamic 

Pressure 

(psf) 

Density 

(Slug/ft 3 ) 

10" 4 

Freq. 

(Hz) 

Reduced 

Freq. 

Mass 

Ratio 

Flutter 

Speed 

Index 

0.600 

669.6 

142.4 

6.35198 

12.74 

0.1132 

171.30 

0.2720 

0.700 

781.2 

124.1 

4.06703 

11.60 

0.0884 

267.54 

0.2540 

0.800 

880.0 

91.1 

2.35279 

9.45 

0.0639 

462.48 

0.2176 

0.850 

935.0 

59.7 

1.36578 

7.63 

0.0486 

796.70 

0.1761 

0.950 

1045.0 

284.42 

5.20904 

15.56 

0.0886 

208.89 

0.3845 

1.05 

1155.0 

358.00 

5.36722 

16.24 

0.0837 

202.73 

0.4313 































































































































































Table 8 CAP-TSD flutter results computed using the nonlinear equation for the extended wing configuration. 


Mach 

No. 

KHSj 

Dynamic 

Pressure 

(psf) 

Density 

(Slug/ft 3 ) 

lcr 4 

Freq. 

(Hz) 

Reduced 

Freq. 

Mass 

Ratio 

Flutter 

Speed 

Index 

0.700 

781.2 

126.0 

4.12930 

11.62 

0.0885 

263.51 

0.2559 

0.800 

880.0 

98.1 

2.53357 

9.93 

0.0671 

429.5 

0.2258 

0.850 

935.0 

63.9 

1.46187 

7.85 

0.0500 

744.33 

0.1822 


Table 9 CAP-TSD flutter results computed using the nonlinear 
equation by varying U for the wing/fuselage configuration. 


Mach 

No. 


Dynamic 

Pressure 

(psf) 

Density 

(Slug/ft 3 ) 

1(T 4 

Freq. 

(Hz) 

Reduced 

Freq. 

Mass 

Ratio 

Rutter 

Speed 

Index 

0.800 

771.3 

93.0 

3.12140 

9.60 

0.0741 

348.60 

0.2196 

0.900 

559.5 

34.2 

2.18364 

6.21 

0.0660 

498.30 

0.1333 


Tab 

equatio 

le 10 CAP-TSD flutter results computed using the nonlinear 
n for different mass ratios with the wing/fuselage configuration. 


Mach 

No. 

Velocity 

(ft/sec) 

Dynamic 

Pressure 

(psf) 

Density 

(Slug/ft 3 ) 

nr 4 

Freq. 

(Hz) 

Reduced 

Freq. 

Mass 

Ratio 

Flutter 

Speed 

Index 

0.900 

646.0 

113.54 

5.44056 

28.48 

0.2623 

200.00 

0.2429 

0.900 

559.5 

34.2 

2.18364 

6.21 

0.0660 

498.30 

0.1333 


[able 11 CF1 

J3D-Euler flutter results computed for the extended wii 

lg configurati 

on. 

Mach 

No. 

Velocity 

(ft/sec) 

Dynamic 

Pressure 

(psf) 

Density 

(Slug/ft 3 ) 

lO" 4 

Freq. 

(Hz) 

Reduced 

Freq. 

Mass 

Ratio 

Rutter 

Speed 

Index 

0.700 

781.2 

118.1 

3.87040 

11.81 

0.0900 

281.14 

0.2477 

0.850 

935.0 

69.1 

1.58083 

8.62 

0.0549 

688.3 

0.1895 

0.900 

990.0 

32.83 

0.66993 

6.24 

0.0375 

1624.2 

0.1306 

1.050 

1155.0 

266.2 

3.99020 

20.77 

0.1070 

272.7 

0.3719 
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Table 12 CFL3D-Euler flutter results computed for 
different mass ratios with the extended wing configuration. 


Mach 

No. 

Velocity 

(ft/sec) 

Dynamic 

Pressure 

(psf) 

Density 

(Slug/ft 3 ) 

10“ 4 

Freq. 

(Hz) 

Reduced 

Freq. 

Mass 

Ratio 

Flutter 

Speed 

Index 

0.850 

760.0 

72.2 

2.50185 

8.81 

0.0690 

434.9 

0.1937 

0.900 

667.0 

48.6 

2.18364 

7.07 

0.0632 

498.30 

0.1583 


Table 13 CFL3D-Navier-Stokes flutter results computed with the extended wing configuration. 


Mach 

No. 

Velocity 

(ft/sec) 

Dynamic 

Pressure 

(psf) 

Density 

(Slug/ft 3 ) 

10" 4 

Freq. 

(Hz) 

Reduced 

Freq. 

Mass 

Ratio 

Flutter 

Speed 

Index 

0.850 

935.0 

80.4 

1.83866 

8.80 

0.0560 

591.8 

0.2044 

0.900 

990.0 

79.8 

1.62881 

8.16 

0.0491 

668.0 

0.2037 

0.920 

1012.0 

111.7 

2.18133 

9.10 

0.0535 

498.8 

0.2410 
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Figure 6 Symmetry plane and wing surface grid used for CFL3D Navier-Stokes calculations. 
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Figure 7 Effects of time step on the 
first generalized displacement computed 
with CAP-TSD using the nonlinear 
equation at Moo=0.6 and q co = 136.8 psf. 


Figure 8 Effects of time step on the 
first generalized displacement computed 
with CAP-TSD using the nonlinear 
equation at Moo =0*9 and qoo =57.6 psf. 


— 100x50x80 



Figure 9 Effects of grid density on the first generalized displacement computed 
using CAP-TSD with the nonlinear equation at Moo =0.9 and qoo =2 1.6 psf. 
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Figure 12 Effect of generalized velocity initial Figure 13 Effect of generalized velocity initial 

conditions on the first generalized displacement conditions on the first generalized displacement 

computed with CAP-TSD using the nonlinear computed with CAP-TSD using the nonlinear 

equation at Moo =0.85 and qoo=64.8 psf. equation at Moo =0.90 and qoo=21.6 psf. 



x/c 


x/c 


Figure 14 Comparison of CAP-TSD and CFL3D- Euler steady pressure 
coefficients for the extended rigid wing configuration at Moo =0.7 and a = 0. 
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Wing Root 


and q o =93.6 psf 



Figure 18 Effects of freestream Mach number (and dynamic pressure) on steady 
pressure coefficients computed using CAP-TSD for the wing/fuselage configuration 
under static aeroelastic deformation with dynamic pressure near the flutter value. 
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(d) Moo =0-95 and qoo=259.2 psf. 
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Figure 18 Effects of fuselage modeling on steady pressure 
coefficients computed using CAP-TSD at M=0.85 and q=64.8 psf. 


31 






y/s = 0.95 







Figure 21 Comparison of sectional lift coefficients computed using CAP-TSD and CFL3D-Euler 
for static aeroelastic cases at (a) Moo =0.7, qoo=115.2 psf and (b) M oo=0-9, qoo=36,0 psf. 




Figure 22 Aerodynamic center location in relation to the first torsion node line computed using CAP-TSD 
for static aeroelastic cases at (a) Moo =0.7, q oo=115.2 psf and (b) Moo=0.9. q«,=36.0 psf. 




Figure 23 CAP-TSD flutter points computed using the linear and nonlinear 
equations compared with Experimental data for the wing/fuselage configuration. 



Figure 24 Effects of fuselage modeling on the predicted flutter boundary using the nonlinear CAP-TSD equation. 



Figure 25 Effects of mass ratio on the predicted flutter boundary using the nonlinear CAP-TSD equation for the wing/fuselage configura 
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Figure 27 Comparison of CFL3D-Euler and CAP-TSD computed flutter 
solutions with experimental data for the extended wing configuration. 
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